REPORT  DOCUMENTATION  PAGE 


AFRL-SR-BL-TR-99- 


Public  reporting  burden  for  this  collection  of  information  is  estimated  to  average  1  hour  per  response 
gathering  and  maintaining  the  data  needed,  and  completing  and  reviewing  the  collection  of  informatii 
collection  of  information,  including  suggestions  for  reducing  this  burden,  to  Washington  Headquarter! 
Davis  Highway,  Suite  1204,  Arlington,  VA  22202-4302,  and  to  the  Office  of  Management  and  Budgr 


1.  AGENCY  USE  ONLY  ( Leave  blank) 


REPORT  DATE 


l w&fo 


lata  sources, 
spect  of  this 
1  5  Jefferson 
503. 


1/29/99 


4.  TITLE  AND  SUBTITLE 

Experimental  and  Numerical  Investigations  on  Failure  Mechanisms  in 
Fiber-Reinforced  Ceramics 


Final  Technical  Report  8-1-96  to  6-30-98 


5.  FUNDING  NUMBERS 

F49620-95- 1-0262 


6.  AUTHOR(S) 
H.  L.  Schreyer 
M.  L.  Wang 


61102F 

2302/BS 


7.  PERFORMING  ORGANIZATION  NAME(S)  AND  ADDRESS(ES) 

H.  L.  Schreyer,  University  of  New  Mexico,  Department  of  Mechanical  Engineering 
Albuquerque,  NM  87131 


8.  PERFORMING  ORGANIZATION 
REPORT  NUMBER 


M.  L.  Wang,  University  of  Illinois  at  Chicago,  Department  of  Civil  &  Materials 
Engineering,  842  West  Tavlor.  Chicago.  IL  60607-7023 _ 


9.  SPONSORING/MONITORING  AGENCY  NAME(S)  AND  ADDRESS(ES) 

Air  Force  Office  of  Scientific  Research/NA 
801  N.  Randolph  Street,  Rm  732 
Arlington,  VA  22203-1977 


10.  SPONSORING/MONITORING 
AGENCY  REPORT  NUMBER 

F49620-95- 1-0262 


12a.  DISTRIBUTION  AVAILABILITY  STATEMENT 

Approved  for  public  release;  distribution  unlimited. 


12b.  DISTRIBUTION  CODE 


13.  ABSTRACT  (Maximum  200  words) 

The  primary  focus  on  this  research  was  to  model  the  evolution  of  microcracks  at  the  interfaces  of  grains  and  inclusions  in 
quasibrittle  materials.  Several  unique  aspects  have  been  developed  as  part  of  the  research.  These  include  the  application  of 
the  material  point  method  to  a  new  class  of  problems  and  the  formulation  of  a  thermodynamically-based  decohesion 
constitutive  equation.  In  addition,  numerical  algorithms  were  formulated  to  handle  three-dimensional  problems,  to  include 
the  effects  of  randomness  of  grain  sizes  and  of  the  strength  of  cement  paste,  and  to  take  into  account  boundary  constraints. 
Each  of  these  aspects  represent  new  or  extensive  improvements  to  the  approaches  available  at  the  time  the  research  was 
initiated. 


17.  SECURITY  CLASSIFICATION 
OF  REPORT 

Unclassified 


18.  SECURITY  CLASSIFICATION 
OF  THIS  PAGE 

Unclassified 


15.  NUMBER  OF  PAGES 

48 


16.  PRICE  CODE 


19.  SECURITY  CLASSIFICATION  20.  LIMITATION  OF 
OF  ABSTRACT  ABSTRACT 


Unclassified 


Standard  Form  298  (Rev.  2-89}  (EG) 

Prescribed  by  ANSI  Std.  239.18 

Designed  using  Perform  Pro,  WHS/DIOR,  Oct  94 


Final  Report  for  the  Period:  1  Aug.  1996  -  30  June  1998 


EXPERIMENTAL  AND  NUMERICAL  INVESTIGATIONS  ON  FAILURE 
MECHANISMS  IN  FIBER-REINFORCED  CERAMICS 


by 

H.L.  Schreyer 

Department  of  Mechanical  Engineering 
University  of  New  Mexico 
Albuquerque,  NM  87111 

and 

M.L.  Wang 

Department  of  Civil  and  Materials  Engineering 
University  of  Illinois  at  Chicago 
842  West  Taylor  Street 
Chicago,  IL  60607-7023 


AFOSR  Grant  F49620-95-1-0262 

(Started  1  April  1995) 


Submitted  to: 

Dr.  Ozden  O.  Ochoa 

Mechanics  and  Materials  Program 
AFOSR/NA 
110  Duncan  Ave.,  B115 
Bolling  AFB,  DC  20332-8080 


January  1999 


1  9990302  0  4  9 


EXECUTIVE  SUMMARY 

The  addition  of  spherical  inclusions  and  fibers  into  a  traditional  ceramic  matrix 
can  provide  enhanced  ductility  and  toughness  over  that  of  the  monolithic  ceramic.  The 
primary  objective  of  the  research  was  to  provide  an  improved  microstructural 
understanding  of  a  fiber-reinforced  composite  material.  To  properly  identify  the 
fundamental  issues,  a  combination  of  expertise  in  experimental  methodology  and  of  the 
formulation  of  advanced  constitutive  equations  is  required.  Experimental  data  shows  the 
dominant  failure  mechanisms  and  the  theoretical  investigation  focuses  on  these  features. 

Preliminary  experimental  data  for  a  particular  alumina  showed  that  failure 
occurred  almost  exclusively  along  grain  boundaries.  The  obvious  conclusion  is  that  the 
properties  of  material  that  constitute  grain  boundaries  must  have  a  dominating  influence 
on  the  failure  mode  and  strength  of  the  ceramic.  A  second  observation  was  that  there  did 
not  appear  to  be  a  microcrack-cloud  zone,  a  result  that  is  contrary  to  some  reports  in  the 
literature.  Instead,  microcrack  branches  form  along  grain  boundaries  with  one 
microcrack  eventually  becoming  dominant  to  the  point  where  all  contact  is  lost  at  which 
time  the  adjacent  microcracks  become  inactive.  The  process  can  be  thought  of  as 
decohesion  in  which  traction  carrying  capability  is  gradually  reduced  to  zero  as  a 
microcrack  evolves  to  a  macrocrack.  This  procedure  is  also  called  material  failure. 

Microcracks  actually  represent  a  strong  discontinuity  in  displacement,  a  feature 
that  is  particularly  difficult  to  represent  both  theoretically  and  numerically.  Since  failure 
along  the  boundaries  of  inclusions  and  fibers  is  also  one  of  an  evolving  microcrack,  the 
key  component  to  being  able  to  analyze  a  whole  class  of  problems  is  that  of  modeling 
decohesion.  Simultaneously,  there  is  the  problem  of  obtaining  a  robust  numerical 
solution  procedure. 

A  team  at  the  University  of  New  Mexico  has  been  in  the  forefront  of  the 
development  of  a  numerical  scheme  applicable  to  a  wide  range  of  problems  in  continuum 
mechanics  called  the  Material  Point  Method  (MPM).  Part  of  this  research  involved 
investigating  the  suitability  of  the  method  for  handling  material  failure.  Simultaneously, 
a  thermodynamically-based  derivation  of  a  decohesion  model  was  developed  to  the  point 
where  the  evolution  of  decohesion  could  be  handled  through  a  constitutive  equation. 
This  precludes  the  need  for  special  elements,  realignment  of  elements  or  of  continuous 
remeshing  as  a  crack  propagates.  Because  a  length  scale  is  a  natural  part  of  a  decohesive 
constitutive  equation,  convergence  with  mesh  refinement  can  be  expected  and  is,  indeed, 
shown.  Schreyer,  Sulsky  and  Zhou  have  provided  a  study  of  ceramic  grains  showing  the 
potential  effect  of  how  damage  may  evolve  with  cooling. 

A  true  representation  of  the  behavior  of  quasibrittle  materials  must  take  into 
account  both  the  three-dimensional  nature  of  the  problem  and  the  randomness  of  grain 
sizes  and  the  strength  of  cement  paste.  Wang  and  Chen  have  provided  both  numerical 
and  experimental  data  to  show  that  such  a  fundamental  approach  can  provide  the 
observed  features  of  failure  including  the  important  effects  of  boundary  constraints. 

It  was  quickly  realized  that  a  comprehensive  study  of  the  effects  of  inclusions  and 
fibers  could  only  be  performed  provided  a  complete  description  of  the  matrix  material  is 
understood.  Therefore,  the  direction  of  the  research  became  more  fundamental  and 
concentrated  on  grains  and  inclusions,  and  their  boundaries.  Results  of  these  efforts  are 
provided  by  three  papers  reproduced  here  as  sections  of  the  final  report. 
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MODELING  MATERIAL  FAILURE  AS 
A  STRONG  DISCONTINUITY  WITH 
THE  MATERIAL  POINT  METHOD 


Howard  L.  Schreyer1,  Deborah  L.  Sulsky2 
and  S.-J.  Zhou1 

*Dept,  of  Mechanical  Engineering 
2 Dept,  of  Mathematics  and  Statistics 
University  of  New  Mexico 
Albuquerque ,  NM  USA  87131 


ABSTRACT: 

A  discrete  constitutive  equation  for  modeling  material  failure  as  a  decohesion 
or  separation  of  material  to  form  two  free  surfaces  is  a  relatively  simple  approach. 
However ,  numerical  simulations  based  on  such  a  model  involve  considerable 
complexity  including  remeshing  if  the  finite  element  approach  is  used.  Here,  a 
basic  formulation  involving  the  simultaneous  application  of  the  continuum  and 
decohesion  constitutive  equations  is  described  together  with  a  numerical  approach 
based  on  the  material  point  method.  Preliminary  results  indicate  that  failure 
propagation  can  be  predicted  at  an  arbitrary  angle  without  the  dispersion  of  the 
crack  front  that  is  often  observed  with  conventional  finite  elements. 

KEY  WORDS:  Decohesion,  localized  deformation,  softening,  material  point 
method ,  strong  discontinuity. 
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1.  Introduction 


The  propagation  of  cracks  through  concrete  is  just  a  manifestation  of  material 
failure.  Numerous  models  of  failure  have  been  proposed  together  with  numerical 
procedures  for  obtaining  solutions  to  the  governing  boundary  value  problem. 
Unfortunately,  the  material  model  for  failure  is  often  entwined  with  the  numerical 
approach  so  that  it  is  difficult  to  determine  which  aspect  is  the  limiting  component 
if  predictions  do  not  match  experimental  data.  Here,  we  attempt  to  carefully 
differentiate  these  two  essential  components  by  first  concentrating  on  a  basic  failure 
model,  and  then  proposing  the  use  of  a  relatively  new  numerical  procedure,  the 
material  point  method.  Preliminary  results  indicate  that  the  difficulties  associated 
with  the  finite  element  method  are  not  present. 

There  are  many  criteria  for  material  failure  but  the  definitions  of  failure  are  often 
vague  or  defined  implicitly  through  each  criterion.  We  take  material  failure  to  mean 
the  process  by  which  two  new  free  surfaces  are  formed,  with  brittle  fracture  as  an 
obvious  example.  However,  there  are  other  forms  of  material  failure  as  exemplified 
by  ductile  rupture,  delamination,  the  breaking  of  grain  boundaries  and  the  pullout  of 
reinforcing  rods  or  fibers.  Our  interest  is  to  represent  all  of  these  phenomena  with  a 
single  model  that  incorporates  the  essential  features  of  the  state  of  stress  or  strain  at 
which  failure  initiates  and  predicts  the  correct  energy  dissipated.  Our  focus  is  not  on 
replicating  the  details  of  failure,  although  this  can  be  done  in  some  cases,  but  on 
predicting  the  effect  of  failure  on  the  far-field  stress  distribution  and  on  structural 
response  as  reflected,  for  example,  by  a  force-deflection  curve.  The  proposed 
approach  is  a  constitutive  equation  that  describes  decohesion.  When  used  with  the 
material  point  method,  which  is  a  relatively  new  computational  method  that  is 
particularly  robust  for  problems  with  large  deformations,  the  proposed  approach  has 
a  simple  structure  in  that  the  decohesion  comes  into  the  analysis  through  the 
constitutive  equation  only.  There  is  no  attempt  to  enforce  the  geometrical 
continuity  of  a  crack.  Instead,  compatibility  is  enforced  in  an  averaged  sense. 

Failure  modeling  involves  both  theoretical  formulations  of  constitutive  equations 
and  numerical  simulations,  and  the  two  aspects  should  be  carefully  delineated. 
However,  the  finite  element  method  has  become  the  method  of  choice  for  the 
majority  of  engineering  applications  so  that  the  formulation  of  the  constitutive 
equations  is  often  tailored  for  use  by  finite  elements;  conversely  limitations  imposed 
by  the  finite  element  method  are  often  interpreted  unjustly  as  a  limitation  of  the 
theoretical  approach.  In  the  following  brief  survey,  we  attempt  to  keep  the 
discussion  of  the  two  phases  distinct  if  at  all  possible. 

A  large  number  of  papers  related  to  failure  have  been  based  on  a  zone  of  softening 
with  an  assumed  width  in  which  a  continuum  constitutive  equation  continues  to  be 
used  [BAZ  84,  BOR  87,  ROT  87,  OLI  89,  DAH  90,  WIE  98].  A  theoretical 
difficulty  with  such  an  approach  is  the  possible  loss  of  ellipticity  and  material 
stability  within  the  band.  When  used  with  finite  elements,  the  band  width  is 
associated  with  the  size  of  the  elements  and  the  accuracy  is  then  limited  when  the 
elements  become  highly  deformed. 
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An  alternative  (discrete)  approach  is  to  consider  material  failure  as  a  strong 
discontinuity  in  displacement  with  traction  related  to  the  discontinuity.  There  is  a 
long  history  in  which  discrete  constitutive  equations  are  postulated  directly  as 
reflected  by  Barenblatt  [BAR  59],  and  Hillerborg  et  al.  [HIL  76].  Feenstra  et  al. 
[FEE  1991]  and  Corigliano  [COR1993]  provide  a  nice  summary  of  previous  models 
and  describe  numerical  methods  based  on  the  use  of  interface  elements.  The  use  of 
discrete  constitutive  equations  has  not  met  with  complete  favor  partially  because 
strong  discontinuities  are  difficult  to  handle  numerically  and  convergence  with  mesh 
refinement  and  mesh  insensitivity  is  difficult  to  show.  The  use  of  interface  elements 
may  require  frequent  remeshing  if  the  crack  surface  propagates  in  a  curved  manner, 
and  double  nodes  which  separate  with  the  evolution  of  decohesion  [SCH,  1992], 
However,  Dvorkin  et  al.  [DVO  90]  provide  a  nice  approach  that  overcomes  many  of 
these  objections  by  handling  discontinuities  at  the  element  level  rather  than 
enforcing  discontinuities  to  be  along  element  boundaries. 

A  fundamentally  different  approach  is  described  in  more  recent  work  by  Simo  et 
al.  [SIM  93]  in  which  the  continuum  constitutive  equation  is  extended  beyond  the 
loss  of  ellipticity  condition  into  the  softening  regime.  They  argue  that  this 
extension  should  be  accompanied  by  distribution  theory  which,  in  effect,  leads  to  a 
strong  discontinuity.  The  theory  has  since  been  extensively  developed  by  Oliver 
[OLI  89],  Simo  and  Oliver  [SIM  94],  Armero  and  Garikipati  [ARM  95],  Larsson  and 
Runesson  [LAR  96],  Oliver  [OLI  96]  and  Armero  [ARM  97].  The  final  result  is 
discrete  constitutive  equations  relating  stress  to  the  discontinuity  in  displacement, 
and  here  also  the  discontinuity  is  handled  at  the  element  or  constitutive  level. 

We  have  opted  for  a  particular  combination  of  these  ideas  in  an  attempt  to 
provide  an  approach  that  is  as  simple  and  as  straightforward  as  possible.  First,  we 
propose  the  direct  introduction  of  discrete  constitutive  equations  with  the  thought 
that  they  should  be  introduced  when  ellipticity  is  lost,  although  a  direct  failure 
initiation  criterion  can  be  used.  No  attempt  is  made  to  model  the  post-crack 
frictional  effects  that  may  occur  with  surfaces  with  rough  cracks  [FEE  91]  although 
such  features  can  be  added.  Second,  the  discontinuity  is  considered  to  be  part  of  the 
constitutive  equation  and  is  applied  in  a  manner  analogous  to  that  of  Dvorkin  [DVO 
91],  Oliver  [OLI  96a,  OLI  96b]  and  Armero  [ARM  97].  A  point  that  is  undergoing 
failure  is  also  considered  to  be  a  material  point  in  the  continuum  so  that  the 
decohesion  and  continuum  constitutive  equations  must  be  simultaneously  satisfied 
subject  to  the  restriction  of  traction  equilibrium.  Third,  we  invoke  the  constitutive 
equation  in  the  material  point  method.  The  arguments  for  the  direct  calculation  of 
the  strong  discontinuity  in  displacement,  which  we  also  call  decohesion,  and  the  use 
of  the  material  point  method  are  summarized  as  follows: 

(i)  We  retain  the  conceptual  simplicity  inherent  with  the  discrete  constitutive 
approach  that  material  failure  does  not  happen  abruptly  but  occurs  smoothly  with  a 
gradual  reduction  in  traction  as  the  displacement  discontinuity  increases. 

(ii)  We  believe  it  is  extremely  difficult  to  evaluate  properties  of  any  constitutive 
equation  in  the  failure  regime.  However,  it  is  probably  easier  to  select  material 
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parameters  for  a  discrete  constitutive  equation  than  for  a  continuum  model  extended 
into  the  softening  regime. 

(iii)  The  discrete  equation  can  be  applied,  if  desired,  at  the  instance  ellipticity  is  lost 
so  that  there  is  a  high  probability  that  well  posedness  can  be  retained  although  a 
stability  analysis  must  be  performed  [SUO  92]. 

(iv)  The  essential  aspects  of  prescribed  stress  at  the  initiation  of  failure  and  prescribed 
energy  dissipation  at  the  end  of  failure  are  automatically  included  in  this  model. 

(v)  Once  decohesion  is  initiated  on  a  surface  of  discontinuity,  the  adjacent  continuum 
tends  to  unload  into  the  elastic  regime,  so  the  computational  simplicity  of  only 
needing  to  combine  decohesion  with  elasticity  covers  the  vast  majority  of  practical 
cases. 

(vi)  The  decohesion  constitutive  equation  can  be  developed  in  a  thermodynamical 
setting,  in  concert  with  many  current  continuum  models,  and  can  include  plasticity, 
damage,  viscoelastic  and  viscoplastic  features  that  are  associated  strictly  with  the 
decohesion. 

(vii)  The  application  of  decohesion  constitutive  equations  in  the  material  point 
method  retains  the  simplicity  of  current  applications  of  strong  discontinuities  at  the 
element  level  in  the  finite  element  method.  However,  double  nodes  or  interface 
elements  are  not  needed  and  there  is  the  additional  potential  advantage  that  mesh 
orientation  and  mesh  distortion  are  not  factors  that  need  to  be  considered. 

(viii)  Following  the  method  outlined  by  Allix  and  Corigliano  [ALL  95,  ALL  96] 
there  is  the  potential  of  relating  the  decohesion  constitutive  equation  to  mixed-mode 
fracture. 

(ix)  Finally,  the  use  of  a  discrete  constitutive  equation  may  still  be  a  suitable  model 
for  diffuse  failure  if  the  primary  objective  is  to  obtain  an  efficient  solution  for  the 
region  away  from  the  failure  zone. 

The  next  section  provides  only  a  brief  description  of  the  material  point  method 
since  the  method  has  been  fully  described  in  previous  papers.  Section  3  describes  the 
basic  structure  of  the  decohesion  model  used  in  our  analysis.  Analytical  and 
numerical  solutions  to  model  problems  [ZHO  98]  including  a  convergence  study  are 
given  in  Section  4  which  is  then  followed  by  conclusions  concerning  the  general 
applicability  of  the  method  for  material  failure  in  general  including  delamination. 


2.  The  Material  Point  Method 

The  material  point  method  [SUL  94,  SUL  95,  SUL  96]  discretizes  a  solid  body 
by  marking  a  set  of  material  points  in  the  original  configuration  that  are  tracked 
throughout  the  deformation  process.  Let  xj,  p  =  l,...,Np  denote  the  current  position 

of  material  point  p  at  time  tn,  n  =  0,  1,  2,  ....  These  material  points  provide  a 
Lagrangian  description  of  the  solid  body  that  is  not  subject  to  mesh  tangling.  Each 
point  at  time  t°  has  an  associated  mass,  mp,  density,  pj,  velocity,  v“,  Cauchy  stress 

tensor,  c“,  strain,  ej,  and  any  other  internal  variables  necessary  for  the  constitutive 
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model.  If  temperature  changes  are  important,  internal  energy  or  temperature  may 
also  be  ascribed  to  the  material  points.  The  material  point  mass  is  constant  in  time, 
insuring  that  the  continuity  equation  is  satisfied.  Other  variables  must  be  updated 
with  reference  to  conservation  of  momentum,  conservation  of  energy,  or  from  the 
constitutive  model. 

To  make  the  computations  tractable,  at  each  timestep  of  a  dynamic  algorithm, 
information  from  the  material  points  is  interpolated  to  a  background  computational 
mesh.  This  mesh  covers  the  computational  domain  and  is  chosen  for  computational 
convenience.  A  particularly  simple  choice  is  a  regular  rectangular  grid.  After 
information  is  interpolated  to  the  grid,  equations  of  motion  are  solved  on  this  mesh 
which  is  considered  to  be  an  updated  Lagrangian  frame.  For  example,  to  solve  the 
momentum  equation  on  the  grid  using  an  explicit  FE  algorithm,  one  must  know  the 
value  of  the  momentum  at  the  beginning  of  the  timestep  at  the  nodal  positions.  The 
nodal  momentum,  m“v“,  is  the  product  of  the  nodal  mass  and  nodal  velocity,  and 
each  is  determined  by  interpolation, 


Np 


p-1 


In  the  above,  Nj(x)  is  the  nodal  basis  function  associated  with  node  I.  In  this  paper, 
Nj(x)  are  the  tensor  products  of  piecewise  linear  functions.  The  internal  forces  are 
determined  from  the  particle  stresses  according  to 


The  quantity  G“p  is  the  gradient  of  the  nodal  basis  function  evaluated  at  the  material 
point  position,  G“p  =  VNi(x)|x0.  The  momentum  equation  is  solved  with  the  nodes 


considered  to  be  moving  with  the  deformation  to  give  nodal  velocities,  v,L ,  at  the 
end  of  this  Lagrangian  timestep  of  size  At, 


m“ 


vf  -  v? 


At 


=  f“ 


(3) 


At  the  end  of  this  Lagrangian  step,  the  new  nodal  values  of  velocity  are  used  to 
update  the  material  points.  The  material  points  move  along  with  the  nodes 
according  to  the  solution  given  throughout  the  elements  by  the  nodal  basis  functions 


(4) 


Similarly,  the  material  point  velocity  is  updated  via 

vj+l  -  vj  +^(v}-  -vf  JN,(xJ).  (5) 

i«l 

The  sums  in  these  last  two  equations  extend  from  1  to  Nn  where  Nn  is  the  number  of 
nodes  in  the  computational  mesh. 

A  strain  increment  for  each  material  point  is  determined  using  the  gradient 
of  the  nodal  basis  function, 

Ae;=^^[G^+(G»v^)T].  (6) 

^  i-1 
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This  strain  increment  is  then  used  in  an  appropriate  constitutive  equation  for  the 
material  being  modeled  to  update  the  stress  at  the  material  point.  Any  internal 
variables  necessary  in  the  constitutive  model  can  also  be  assigned  to  the  material 
points  and  transported  along  with  them.  Once  the  material  points  have  been 
completely  updated,  the  computational  mesh  may  be  discarded  and  a  new  mesh 
defined,  if  desired,  and  then  the  next  timestep  is  begun. 

The  material  point  method  has  several  advantages.  The  Lagrangian  description 
provided  by  the  material  points  can  undergo  large  deformations  without  mesh 
tangling.  Since  the  computational  mesh  is  under  user  control,  it  can  be  chosen  so 
that  reasonable  timesteps  may  be  taken  in  this  Lagrangian  frame.  Usually,  the 
timestep  is  restricted  by  the  CFL  condition  for  an  explicit  algorithm,  where  the 
critical  timestep  is  the  ratio  of  the  mesh  size  to  the  wave  speed.  Note  that  this 
condition  depends  on  the  more  favorable  mesh  spacing,  not  the  material  point 
spacing.  Since  equations  are  solved  in  an  updated  Lagrangian  frame  on  the  FE  mesh, 
the  nonlinear  convective  terms  troublesome  in  Eulerian  formulations,  are  not  an 
issue.  Finally,  the  material  points  transport  material  properties  and  internal 
variables  without  error. 


3.  Discrete  Constitutive  Equation  For 
Decohesion 

3.1  The  Theoretical  Model 

We  define  the  initiation  of  material  failure  as  the  time  when  a  material  point 
first  experiences  a  discontinuity  in  displacement  but  continues  to  function  as  a  point 
in  a  solid  continuum.  A  collection  of  such  points  in  a  neighborhood  defines  a 
failure  surface,  T.  Although  the  material  manifestation  is  a  single  surface,  one 
observes  spatially  two  surfaces,  Tv  and  TL,  as  sketched  in  Fig.  1.  Each  dotted  line 
illustrates  points  in  space  identified  with  a  single  material  point  which  can  be 
considered  associated  with  any  one  of  the  spatial  points  on  the  line.  The  sketch 
illustrates  the  material  surface  as  a  thick  line  between  the  two  spatial  surfaces  but  if 
a  Lagrangian  description  is  used,  the  material  surface  may  be  at  a  totally  different 
location.  Failure  is  said  to  be  complete  when  traction  can  no  longer  be  sustained  on 
the  material  surface,  i.e.,  the  spatial  surfaces  no  longer  have  any  ligaments 
connecting  them  even  though  one  point  on  each  surface  is  identified  as  a  single 
material  point.  The  discontinuity  in  displacement  is  called  decohesion. 

Here,  we  present  a  development  of  discrete  constitutive  equations  using 
thermodynamics  as  a  framework  with  the  result  that  the  dissipation  inequality  is 
automatically  satisfied.  The  approach  entails  two  essential  assumptions  consisting 
of  (1)  the  form  of  the  free  energy,  and  (2)  the  form  of  the  evolution  equations.  Each 
assumption  leads  to  a  different  model  which  can  only  be  tested  by  solving  a  problem 
for  which  either  qualitative  or  quantitative  data  exist.  To  allow  for  the  presentation 
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of  different  models  in  a  convenient  manner,  we  present  the  general  framework  first, 
and  then  show  the  implications  inherent  in  specific  assumptions. 


rL 


(a)  Decohesion  in  spatial  configuration  -  two  surfaces. 


n 


(b)  Decohesion  in  material  configuration  -  surface 
with  discontinuities 

Fig.  1.  Material  failure  as  represented  by  two  spatial  surfaces,  rv  and  TL, 
or  one  material  surface,  T. 

The  approach  is  analogous  to  what  one  might  use  for  a  rigid-plastic  continuum 
for  which  the  elastic  part  of  the  response  is  ignored,  i.e.,  the  total  strain  and  the 
plastic  strain  are  identical.  The  internal  strain  energy  does  not  exist  and  the  stress 
must  be  provided  by  the  solution  to  a  boundary  value  problem.  However,  there 
remains  a  contribution  to  the  free  energy  associated  with  hardening  and  evolution 
equations  for  plasticity  variables. 

Consider  a  situation  where  loads  are  applied  to  a  body  which  is  continuous  except 
on  a  material  failure  surface,  T,  which  displays  a  strong  discontinuity,  or  decohesion, 
ud  =  [u].  Any  point  on  the  suface  is  also  a  point  in  the  continuum  which  is 
assumed  to  be  governed  by  linear  elasticity  so  that  the  stress,  a,  and  strain,  e,  are 
linearly  related  by  the  elasticity  tensor,  E: 

a  =  E:e  (7) 

If  n  denotes  the  normal  to  the  surface,  then  the  traction,  T,  is  given  by  T  =  G-n. 
The  rate  at  which  power  is  being  added  to  the  surface  by  this  traction  is  x  •  ud  in 
which  a  superposed  dot  denotes  a  derivative  with  respect  to  time.  We  postulate  that 
the  free  energy  per  unit  surface  area  consists  of  an  initial  energy,  U0,  due  to  residual 
stresses  that  resulted  from  the  curing  process,  and  a  term,  Ud,  which  represents  the 
effect  of  decohesion: 
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U  =  U0  -Ud(u)  (8) 

The  use  of  the  negative  sign  is  meant  to  suggest  that  normally  energy  is  provided  by 
the  original  material  to  the  decohesion  process.  The  parameter,  u,  is  a  scalar 
representation  of  the  state  of  decohesion.  The  specific  choice  for  Ud  (positive)  is  part 
of  a  particular  model  The  decohesion,  up,  is  viewed  as  “permanent”  decohesion  and 
is  introduced  similarly  to  plastic  strain  in  elasto-plastic  continuum  models.  In  the 
absense  of  elasticity,  up  -  ud. 

If  only  up  and  u  are  considered  to  be  the  primary  variables  describing  decohesion, 
the  dissipation  rate  is 

Ds  =  xud -U  =  x-iip+xu  with  T  =  ^=r-  (9) 

The  generalized  traction,  X,  is  conjugate  to  u.  Instead  of  the  traction  starting  at 
zero,  as  it  does  for  some  existing  discrete  models  [NEE  87,  COR  93],  we  visualize 
that  when  the  failure  process  starts,  ud  =  0,  up  =  0,  and  the  traction  is  the  initial 
vector,  T0  which  depends  on  the  path. 

We  parameterize  the  development  of  decohesion  through  a  single,  dimensionless 
monotonically  increasing  variable,  X ,  and  the  evolution  equations 

iip  =  Xmc  u  =  Xme  (10) 

in  which  me  and  me  denote  evolution  functions  that  depend  on  X  and  x.  me  is  a 
vector  whose  inner  product  with  X  is  assumed  to  be  positive,  semi-definite.  If  we 
introduce  an  effective  traction 

xe=xme  (11) 

then  the  dissipation  rate  becomes 

Ds=i[xe+xme]  (12) 

To  ensure  the  dissipation  is  positive,  define  a  decohesion  function  as  follows: 

Fd  =r  +fme  -F0  F0  >  0  (13) 

The  function  has  been  constructed  in  the  usual  manner  so  that  Fd  is  negative  when 
the  tractions  are  zero.  We  assume  decohesion  does  not  occur  unless  Fd  =  0  in  which 
case  the  dissipation  rate  becomes  Ds  =  XF0 ,  a  positive  scalar  (and  Fd  >  0  is  not 
permitted).  The  evolution  equations  can  now  be  interpreted  as  parameterized  in  terms 
of  dissipation  which  is  a  monotonically  increasing  parameter.  The  total  dissipated 
energy  is  simply  D  =  AF0  which  depends  only  on  the  value  of  X  and  not  on  the  path 
followed  to  achieve  a  given  value  of  X. 

Next,  we  consider  the  decohesion  condition,  Fd  =  0.  At  the  initiation  of 
decohesion,  X  =  x0  and  mc  =  mj .  If  we  assume  x  =  0  at  the  initiation,  then 

F0=xs=x0*m5  (14) 

and  the  decohesion  condition  reduces  to 

xc  =  V0  -Tmc  (15) 

Typically,  x  increases  to  the  point  where  X®  goes  to  zero  and  X  =  xs  =  X£  /  mse  which 
is  defined  to  be  separation.  The  values  of  up  and  u  at  separation  are  denoted  by 
usp  and  us ,  respectively.  Unless  there  is  a  load  reversal  which  brings  the  two  spatial 
surfaces  back  into  contact,  it  is  assumed  that  xe  remains  zero  after  separation. 

With  the  use  of  (8),  the  stored  surface  energy  lost  due  to  separation  is 
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*U=-Ud(us) 


(16) 


The  total  energy  per  unit  area  that  must  be  provided  to  cause  total  separation  is 
variously  called  the  fracture  energy,  or  the  energy  of  separation,  <J>F,  and  consists  of 
the  sum  of  the  stored  energy  and  the  dissipated  energy: 

®p=®u+*D  O7) 

Because  is  negative,  the  fracture  energy  is  less  than  the  dissipation. 

In  summary,  the  resulting  set  of  constitutive  equations  in  rate  form  becomes 

(i)  c  =  E:  e  Continuum  elasticity 

(ii)  x  =  a  •  n  Traction  equilibrium 

(iii)  iip  =  Xm*  u  =  Xmt  Evolution  equations 


(iv) 


Constitutive  equation 


(18) 


(v)  Fd  nr  -(xS-xme)  =  0  Consistency 


(vi)xe=x-mc  andT5=x0*m5 

Two  additional  assumptions  remain  to  completely  formulate  constitutive 
equations;  (i)  the  form  of  the  evolution  equations  for  the  evolution  functions  me  and 
me,  and  (ii)  the  form  of  the  function  Ud  which  provides  a  constitutive  relation 
between  x  and  u.  Even  slight  changes  in  the  forms  of  these  assumptions  can  have 
significant  effects  on  predictions.  Since  the  only  possible  way  to  evaluate  the 
suitability  of  decohesive  constitutive  equations  is  indirectly  through  comparisons  of 
solutions  to  problems  with  features  provided  by  experimental  data,  we  consider  (18) 
to  be  the  basic  format  and  provide  different  models  based  on  plausible  assumptions. 
The  results  of  choosing  specific  forms  for  me,  me  and  Ud  are  given  next. 


3.2  Model  1:  Associated  Evolution  Equations 

Here  we  are  more  specific  in  our  formulation  of  the  decohesion  model.  In  the 
theoretical  formulation,  it  is  most  convenient  to  use  dimensional  parameters  so  that 
physical  interpretations  can  be  easily  made;  conversely,  for  numerical 
implementation  of  the  theory,  dimensionless  variables  should  be  used.  With  these 
objectives  in  mind,  we  choose  X  to  be  dimensionless  and  consider  u  to  have  the 
dimension  of  length  (in  analogy  with  the  decohesion  ud).  We  choose  the  evolution 
functions  and  the  decohesion  energy  to  be  of  the  following  forms: 


me  =  un 


A,  -x 


mc 


;u0 


ud=u0 


(u/u,r‘ 


(19) 


0  (t-A,,  -T)1'2  “  ~“°  Wi-W°  (q  +  1) 

in  which  Ad  is  taken  to  be  a  positive  definite  (dimensionless)  tensor  whose 
components  are  material  parameters,  as  is  q  >  0 .  Additional  material  parameters 
(constants)  are  the  reference  decohesion  scalar,  u0,  and  the  reference  surface  energy, 
U0.  We  define  a  reference  scalar  traction,  x0,  by  the  relation  U0  =u0x0.  The 
immediate  result  of  these  choices  is  that 


Te  =  U0(X-Ad-X)u 


u  =  Xun 


x  =  'c0(X)q 


(20) 
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and  the  evolution  functions  are  obtained  as  derivatives  of  the  damage  function  with 
respect  to  the  corresponding  conjugate  variables,  i.e.,  the  evolution  functions  are  said 
to  be  “associated”: 


(21) 


At  this  point  we  consider  a  two-dimensional  formulation  with  n  denoting  the 
normal  to  the  failure  surface  and  t  a  unit  tangent  vector  as  indicated  in  Fig.  1. 
Corresponding  components  of  x  are  Tn  and  xt,  respectively.  If  the  traction  consists 
only  of  a  normal  component,  specify  the  failure  initiation  value  as  x^  and,  similarly, 
let  Ttf  denote  the  failure  initiation  value  for  a  purely  shear  case.  One  approach  for 
incorporating  these  failure  initiation  conditions  is  to  choose  the  components  of 
with  respect  to  this  local  basis  as  follows: 


[AJ  =  t02 


1 

CO2 

0 


(22) 


The  consequences  of  this  choice  are 


D  =XU0 


Te  =  TJ0Te 


Tc  = 


Fd=U0Fd* 

(x-Ad-x)1/2 


f;  =  xe*-(i-x*) 

tt2 


(23) 


where  xe'  is  a  dimensionless  effective  traction,  x  *  is  a  dimensionless  form  of  x , 
and  Fd  is  a  dimensionless  decohesion  function.  In  deriving  these  equations,  the 


identities  F0  =  x£  =  U0  have  been  used.  The  identities  follow  from  the  conditions  at 


the  initiation  of  decohesion  that  u  =  0  and  xe*  =  1 . 

We  note  that  the  decohesion  condition  (Fd  =  0)  reduces  to  Tc*  =  1  -  x  * .  As 


decohesion  occurs  u  increases  and  xc*  decreases  to  zero  when  u  =  u0 .  Therefore  u0 
can  be  interpreted  as  the  value  of  u  at  which  separation  occurs.  In  the  post 
separation  regime,  u  >  u0,  the  decohesion  condition  is  xc*  =  0 .  Finally,  we  give  an 


alternative  form  for  the  decohesion  evolution  function: 
U0/xn 


me 


Te'"^X; 


^-n  +  rH) 

of 


(24) 


Since  the  dissipation  rate  is  A,U0,  the  energy  dissipated  per  unit  surface  area  at 
any  moment  is  simply  A,U0.  From  (20),  X  =  1  at  separation  so  the  maximum 
energy  dissipated  is  U0  which  provides  a  physical  interpretation  and  a  method  for 
determining  this  particular  parameter.  The  formulation  implies  the  dissipated  energy 
is  independent  of  path  which  is  generally  not  representative  of  real  materials.  The 
final  value  of  the  stored  energy  is  obtained  by  substituting  u  =  u0  in  (19)  and  using 
(16).  In  summary,  the  final  dissipated,  stored  and  total  failure  energies  are: 


<I>D=U0 


- _ l _ 

u  (q  +  1) 


Uft 


3 

F  (q  +  1) 


Un 


(25) 


Suppose  a  pure  opening-mode  path  is  followed,  or  xt  =  0.  Then  it  is  easily 
shown  that  un  =  u(x0/xnf).  If  x0  is  chosen  to  equal  ,  then  u  equals  un  for 
Mode  I.  Experimental  data  obtained  from  a  pure  shear  mode  test  can  then  be  used  to 
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assess  the  adequacy  of  the  model  in  a  process  similar  to  that  used  to  evaluate  Mises 
plasticity.  The  limitation  of  a  single  value  of  dissipated  energy  can  be  circumvented 
by  using  a  multifunctioned  decohesion  surface.  The  development  of  such  a  surface  is 
beyond  the  scope  of  what  we  wish  to  achieve.  For  the  given  model,  the  required  data 
are  xnf,  xtf,  U0  and  q.  We  can  choose  x0  =  and  then  u0  =  U0  /  x^  to  provide 
values  for  all  of  the  parameters. 

Sometimes  nonassociated  models  are  required  to  provide  a  better  fit  with 
experimental  data  including  observations  on  the  mode  of  failure.  Next  we  give  an 
example  of  how  a  particular  nonassociated  model  can  be  constructed. 

3.3  Model  2:  Nonassociated  Evolution  Functions 

Suppose  we  retain  all  aspects  of  the  previous  model  with  the  exception  that  the 
evolution  equation  for  the  permanent  decohesion  is  in  the  normal,  or  opening, 
direction  irrespective  of  the  state  of  traction: 

uj  ss  uf  =  0  m*  =  me  n  (26) 

We  retain  the  previous  expression  for  x*  and  the  decohesion  function.  Therefore,  the 
dissipation  rate  for  normal  mode  decohesion  must  be  evaluated  specifically  from  the 
following  equation: 

D,.n  =XUP+TU  =X[t„uS  +  xu0]  (27) 

which  will  be  less  than  that  obtained  with  the  associated  rule  (sometimes  called  the 
Principle  of  Maximum  Dissipation)  if  xt  is  not  zero  for  at  least  part  of  the 
decohesion. 

A  corresponding  development  for  Mode  II  (pure  shear)  evolution  can  be  obtained 
by  merely  replacing  normal  components  of  traction  with  shear  components. 


4.  Numerical  Application 

4.1  Incorporation  with  the  Material  Point  Method 

In  general,  there  is  no  need  to  determine  the  actual  shape  of  the  deformed  material 
element  associated  with  each  material  point.  However,  when  material  separation 
occurs,  there  is  a  need  to  consider  the  effect  on  the  strain  field  over  the  material 
element  (compatibility).  For  small  deformations,  which  would  be  the  case  normally 
for  quasibrittle  materials  and  for  small  rotations,  the  original  configuration  can  be 
used.  Typically,  each  cell  with  the  material  point  method  is  chosen  to  be  a  square 
element  with  each  side  of  length  h  and  the  element  associated  with  each  material 
point  can  also  be  chosen  similarly  to  be  a  square  of  size  hp  =  h  /  ^n^  where  np  is 

the  number  of  material  points  per  cell.  Over  each  material  element,  the  increment  in 
strain,  Ae,  is  assumed  to  be  constant  as  is  the  increment  in  decohesion,  Aud,  as 
indicated  in  Fig.  4  which  also  shows  the  unit  vector  m  =  Aud  /|Aud|.  For  future 
use,  define  the  opening,  Mn,  shear,  Mt,  and  failure,  Mm,  tensor  modes  as  follows: 
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(28) 


M0=n®n  Mt  =y(n<S>t  +  t®  n)  Mm  = -(n®m  +  m  ®  n) 

We  note  that  Mra  reduces  to  Mn  and  Mt  when  m  =  n  and  m  =  t,  respectively. 

For  a  given,  time  increment,  if  the  total  (average)  strain  increment,  Ae,  is 
considered  fixed,  the  result  of  the  decohesion  is  that  the  effective  strain  increment  in 
the  remaining  part  of  the  material  in  the  element  must  be  reduced  (relaxed)  by  what 
might  be  called  a  decohesion  strain  increment,  Aed,  which  satisfies  a  weak  form  of 
the  compatibility  condition 

J  AeddV  =  J  AudMmdA  Aud  =  Audm  (29) 

tip  and 

in  which  dV  and  dA  denote  differentials  of  volume  on  the  material  element,  Qp,  and 
of  area  on  the  decohesion  surface,  3 £2d,  respectively.  The  magnitude  of  the 
decohesion  increment,  which  is  in  the  direction  of  m  by  definition,  is  Aud.  With  the 
assumptions  that  the  decohesion  and  strain  are  constant  over  each  material  element, 
the  result  is  the  following  expression  relating  the  “relaxation”  or  “decohesion”  strain 
increment  to  the  increment  in  decohesion: 

Ae“=^Mm  Le=^  (30) 

The  effective  length,  Le,  is  merely  the  ratio  of  the  element  volume  to  the  area  of  the 
decohesion  surface  within  that  element.  For  the  two-dimensional  case  illustrated  in 


Fig.  2,  the  effective  length  is 

Le  = — 0<a  <4  (31) 

e  cosap  p  4  v  ' 

with  a  corresponding  formula  for  an  angle  measured  with  respect  to  the  other  side  of 
the  material  element  if  Op  >  rc/4. 


Fig.  2.  A  typical  material  element  with  decohesion. 
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4.2  Solution  Algorithm 

The  constitutive  equations  subroutine  is  invoked  with  the  total  strain  increment, 
Ae,  prescribed,  with  the  total  decohesion  equal  to  the  plastic  decohesion,  and  with  n 
given  and  assumed  fixed.  It  is  computationally  more  efficient  to  define  an  alternative 
mode  vector,  m*,  from  which  the  mode  vector,  me,  and  an  alternative  tensor  mode, 
M*,  are  easily  determined.  Let  0^,  and  denote  the  values  of  0,  up  and  A., 

respectively,  at  the  end  of  the  previous  step.  The  requirement  is  to  solve  the 
following  set  of  nonlinear  equations: 

A0*  =  E:  Ae 

At*  =  A0*  *  n 

me  =  u0m* 

*.=<*«  n 
Aup  =  AAjne 

x  =  Xp,  +  At  (32) 

X  —  A.^  +  AX 

M’  =-yTp-(n®m‘  +m'  <8>n) 

L  L' 

F;  =  r--(i-T*)  =  o 

0  =  0“-  A Xam 


Aa  =  Aa"  -  E:  Ae” 
n 


am  =  E:  M‘ 

Ae*1  =  AXM' 

Ax  =  Ax1  -  AXxn 
up  =  up  +  Aup 

x*  =  r 


The  first  step  is  to  assume  that  no  decohesion  occurs,  to  obtain  a  trial  stress  and 
traction: 


0^=0^+  A0*  x"  =  0"  ■  n  (33) 

and  then  determine  the  value  of  the  damage  function,  Fdtt ,  for  this  trial  traction  and 
the  existing  value  A,pr.  If  <  e  the  step  is  purely  elastic  with  no  additional 
decohesion  and  no  further  action  is  required.  If  the  inequality  is  not  satisfied,  the 
decohesion  variables  must  be  incremented. 

Next  we  describe  a  one-step  algorithm  which  enforces  the  requirement  Fd  =  0  to 
order  (AA-)3 .  Perform  a  Taylor  expansion  of  Fd  about  the  trial  state: 

Fd  =  a(AA.)2  +  2bAA-  +  c  +  0(AA,)3  (34) 


in  which  the  last  term  indicates  the  order  of  the  remainder  and 


a“2 


1  3^ 
dX2 


2b  =  ^ 
ZD  dX 


c  =  E? 


(35) 


We  choose  AA-  =  AXy  and  AA.  =  AA^  to  be  the  solutions  to  the  first-order  and  second- 
order  equations,  respectively;  i.e.,  2bAA^  +  c  =  0  and  a(AA.2)2  +  2bAX2  +  c  =  0  or 

AX,=~^-  A X2=±±^ziL  .  (36) 

with  the  sign  chosen  so  that  in  the  limit  of  infinitesimal  AX  we  have  AX,  =  AX,. 
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Consider  the  case  when  the  model  choice  of  (20)  is  used  and  the  Taylor  expansion 
is  applied  to  the  dimensionless  damage  function  Fd .  It  follows  from  (35)  that 

c  =  (x'T-(l-XV) 


2b-^fl+^ 


(37) 


Za_3A.  dx2 


3x  .  3xe’ 
dX  +  dx 


in  which  all  terms  are  to  be  evaluated  at  and  the  trial  value  of  the  traction.  Note 
that  when  q  =  1,  the  last  term  in  the  expression  for  a  is  zero.  With  a  modest  amount 
of  algebra  involving  (32),  it  follows  that 


3x 

dX 


—  — X_ 


3x1 

3t 


m 


and  we  immediately  have  the  relation 


2b  =  - 


m  x, 


-+qM 


(38) 


(39) 


Next,  we  proceed  to  obtain  the  coefficient  a  in  (37).  We  utilize  (32)  and  (38)  to 
obtain 
dH'" 


dx  .<, 
which  leads  to 
3t  d2x‘‘  dx 
dX  dx2  dX 
Finally 

32x  _  3xm  _ 
dX2  dX  ~ 


[(xm  -Ad  -xm)-(xm  -m')2] 


to  xc- 


„  I?.  3M  _  _  i?.  u0  (  ^  3m  \ 

"nE-ir-"nE-L:(,1®'3r) 


in  which  the  minor  symmetry  of  E  has  been  used.  We  note  that 
3m  ‘  _  3m‘  3x  _  1  32x'”  _  1  M  T 

u 


dX 


dx  dX 


(40) 


(41) 


(42) 


(43) 


with  the  result 

3x1  32x  [(tm  •  Aa  -x„)-(xm  m*)2] 

3x  3 A.2  x2  x'- 

a  result  identical  to  that  of  (41).  Therefore 

_  [(xn  -Ad  -xm)-(xm  m')2]  q(q-l),^ 

a_  xTP 

Even  with  this  rather  simple  form,  preliminary  numerical  results  indicate  that  the 
use  of  the  first-order  equation  for  in  (36)  is  sufficiently  accurate  and  the  extra 
computations  required  to  get  a  is  not  justified. 


(44) 


(45) 


4.3  Separation 

The  procedure  outlined  above  holds  until  separation  occurs  as  indicated  by  X  >  1 . 
The  given  algorithm  can  be  used  also  for  separation  with  the  revised  damage  function 
Fd  =  Te*  so  that  Fd  =  0  enforces  the  separation  condition  that  T  =  0  and  At  =  0.  To 
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prevent  the  possibility  of  numerical  problems  when  Fd  is  close  to  zero,  the 
decohesion  criterion  is  applied  for  FJ  >  £  and  the  separation  condition  is  used  when 
Fd*  <  e  where  e  is  a  small  positive  number  (typically  0.01).  The  value  of  the  stress 
component,  att=tat,  is  automatically  adjusted  through  the  equilibrium 
equations. 

4.4  Example  Solution 

The  application  given  in  this  section  is  restricted  to  a  case  for  which  (i)  the 
initiation  of  failure  is  given  direcdy  by  a  failure  criterion  rather  than  by  a 
discontinuous  bifurcation  analysis,  and  (ii)  the  orientations  of  the  surfaces  of 
decohesion  are  known,  a  priori.  The  purpose  of  these  numerical  analyses  is:  (a)  To 
illustrate  that  the  use  of  jump  in  displacement  as  an  internal  variable  together  with  a 
weak  implementation  of  compatibility  provides  a  simple  and  useful  algorithm  in  the 
MPM,  (b)  to  show  that  the  MPM  does  not  exhibit  the  finite  element  pathologies 
associated  with  distorted  meshes  and  instabilities  with  the  result  that  additional 
features  such  as  enhanced  strains  are  not  required,  and  (c)  to  illustrate  by  example  that 
the  material  point  method  does  not  exhibit  the  orientation  effect  often  seen  with 
finite  elements  when  discontinuities  are  allowed  to  propagate  at  various  angles  to  the 
mesh  sides. 

In  this  example  a  material  composed  of  grains  with  isotropic  elastic  properties 
but  anisotropic  thermal  properties  is  cooled  from  room  temperature  (20°  C)  to  -50° 
C,  at  which  point  the  temperature  is  held  fixed.  The  cooling  rate  is  -10°  C/|is.  The 
coefficient  of  thermal  expansion  is  2.25x lO'V’C  along  the  1-1  material  axis  of  each 
grain  but  an  order  of  magnitude  smaller,  2.1xlO*V°C,  along  the  2-2  material  axis. 
Several  grains  of  this  material  are  made  into  a  composite  material  where  the  grains 
are  assembled  into  the  composite  with  random  orientation. 

Figure  3  shows  a  1  mm  square  sample  of  grains  marked  by  material  points.  An 
arrow  indicates  the  1-1  material  direction  of  each  grain.  The  grains  have  isotropic 
elastic  properties,  with  a  Young's  modulus  3.28  GPa,  Poisson's  ratio  0.3,  and 
density  1.9  kg/m3.  The  interstitial  region  has  the  same  material  properties  and  is 
also  represented  by  material  points  (Fig.  4),  however,  these  points  are  allowed  to 
follow  a  non-associated  decohesion  constitutive  model  in  which  the  mode  is  normal 
to  the  grain  boundary.  The  decohesion  constitutive  properties  are  as  follows:  the 
peak  normal  traction  is  Tnf  =  5.5  MPa,  q  =1,  and  U0  =  0.055  J/m2.  The  peak  stress 
in  shear  was  not  required  because  of  the  nonassociated  evolution  equation. 
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Fig.  3.  Sample  of  anisotropic  grains  with  the  direction  of  the  1-1  material 
axis  indicated  with  an  arrow 


The  normal  to  a  grain  boundary  is  computed  by  assigning  a  scalar  color,  with  the 
value  one,  to  the  material  points  making  up  the  grain,  interpolating  that  color  to  the 
background  mesh  and  then  taking  the  gradient.  The  gradient  gives  the  inward  normal 
to  the  grain.  Figure  5  shows  an  example  of  the  computed  normals  for  a  grain.  TTie 
computational  mesh  is  square  with  side  length  of  25  microns. 


Fig.  4.  Interstitial  region  is  also  discretized  using  material  points. 

As  the  material  cools,  the  anisotropic  thermal  properties  cause  the  contraction  to 
be  nonuniform.  The  effect  of  the  variability  in  the  stress  field  is  seen  in  Fig.  6 
where  contours  of  the  jump  in  displacement  are  shown  at  various  times.  Decohesion 
starts  in  the  upper  right  quadrant  of  the  sample,  and  the  location  of  the  initial 
decohesion  has  the  largest  value  throughout  the  computation.  However,  the 
decohesion  does  not  propagate  from  the  maximum  all  the  way  across  the  domain, 
instead  a  different  path  on  the  left  side  of  the  sample  seems  more  likely  to  span  the 
domain.  This  calculation  indicates  the  complex  interaction  of  geometry  and  material 
properties  that  make  predicting  the  path  of  crack  propagation  so  difficult. 
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0  0,2 
Fig.  6 
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.  Grayscale  contour  plot  of  the  jump  in  displacement  at  time 
(a)  2|ls  and  (b)  5 (is. 
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5.  SUMMARY 


A  rigid,  plastically  softening  decohesion  model  has  been  combined  with 
continuum  elasticity  and  traction  continuity  at  a  material  failure  surface  to  provide  a 
relatively  simple  description  for  the  evolution  of  material  failure.  When  incorporated 
in  the  material  point  method,  the  result  is  a  constitutive  equation  subroutine  that  is 
similar  to  softening  plasticity.  A  length  parameter  associated  with  a  material 
element  provides  a  mechanism  for  ensuring  convergence  with  mesh  refinement.  As 
the  tip  of  a  failure  surface  propagates  through  the  mesh,  the  formulation  inherent 
with  the  material  point  method  appears  to  preclude  the  diffusion  of  the  crack  tip,  a 
feature  often  seen  with  conventional  finite  elements.  Further  investigations 
involving  the  propagation  of  curved  cracks  is  necessary  to  determine  whether  or  not 
the  proposed  method  is  general  and  robust.  Nevertheless,  in  light  of  the  long  history 
of  complex  numerical  analysis  in  connection  with  crack  propagation,  we  believe  the 
simplicity  of  the  decohesion  formulation  in  the  material  point  method  holds 
considerable  promise  for  development  into  a  general  method  for  predicting  material 
failure. 
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Abstract 

Experimental  and  numerical  simulations  are  presented  in  order  to  define  the 
failure  mechanism  of  cementitious  materials. 

Triaxial  tests  on  95  mortar  cylinders  show  there  is  a  change  in  the  failure 
plane  orientation  as  confining  pressure  increases.  The  failure  is  characterized 
by  a  distinct  failure  surface  or  surfaces  that  change  orientation  from  nearly 
vertical  to  angled  as  confinement  increases.  Cap  blocks  are  used  to  relieve  the 
friction  between  the  platens  and  the  specimens  under  triaxial  loading.  The 
same  conclusion  is  obtained  by  numerical  analysis  that  includes  a  realistic  3- 
dimensional  meso-structure  simulation.  A  modified  Gauss  point  method  is 
applied  to  guarantee  the  efficiency  and  accuracy  of  the  finite  element 
analyses.  The  Drucker-Prager  criterion  and  the  low  tensile  strength  criterion 
are  used  to  simulate  the  failure  mechanism. 

Key  words:  Triaxial  test,  random  aggregate  distribution,  modified  Gauss 
point,  failure  mechanisms 


1  Introduction 

Triaxial  tests  have  shown  that  the  angle  of  the  failure  plane  for  geologic 
materials  changes  with  increasing  confining  pressure  (Wawersik,  1971; 
Murrell,  1965).  However,  no  detailed  data  for  concrete  has  shown  that  there  is 
a  similar  relationship  between  the  confinement  and  the  change  in  the  angle  of 
the  principle  failure  plane.  Palaniswamy  (1972),  commenting  on  triaxial  tests 
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on  cement  paste,  mortar,  and  concrete,  stated  that  two  distinct  modes  of 
behavior  were  observed,  specifically,  a  splitting  failure  at  low  confining 
pressures,  and  a  crushing  failure  at  large  confining  pressures.  Zimmerman 
(1977)  found  that  the  failure  plane  is  always  parallel  to  the  direction  of 
maximum  compressive  stress.  The  apparatus  used  in  his  tests  employed 
platens  to  apply  lateral  pressure  to  the  specimens.  These  platens  may  have 
constrained  the  true  failure  modes  of  the  specimens.  Tests  by  Gerstle  (1980), 
using  several  different  loading  devices,  showed  that  boundary  constraints 
from  the  use  of  platens  for  lateral  loading  inhibit  transverse  deformations. 

In  this  paper,  both  experimental  and  numerical  studies  show  changes  in  the 
orientation  of  the  failure  plane  as  confining  pressure  increases. 

Experimentally,  it  is  essential  to  relieve  the  friction  between  the  platens 
since  this  friction  severely  affects  the  development  of  the  failure  mechanism 
and,  consequently,  its  orientation.  Moreover  the  failure  is  characterized  by  a 
distinct  failure  surface  or  surfaces  that  change  orientation  from  nearly  vertical 
to  angled  as  confinement  increases.  The  same  conclusion  is  obtained  by 
numerical  analysis  that  includes  a  realistic  3-dimensional  meso-structure 
simulation.  A  modified  Gauss  point  method  is  applied  to  guarantee  the 
efficiency  and  accuracy  of  the  huge-scale  finite  element  analyses  (148,500 
lattices).  Based  on  comparison  of  the  experimental  and  numerical  results,  the 
strength  ratio  between  cement  paste  (including  interface)  and  coarse  aggregate 
is  randomly  distributed  between  0. 1  to  0.4.  The  Drucker-Prager  criterion  and 
the  low  tensile  strength  criterion  are  used  to  simulate  the  failure  mechanism. 


2  Experimental  Technique 

Triaxial  tests  were  performed  with  confining  pressures  ranging  from  0  to  56 
Mpa.  A  total  of  95  cylinders  were  tested  in  triaxial  compression.  Evaluating 
the  effect  of  the  methods  used  to  relieve  friction  between  the  load  platens  and 
the  specimens  was  essential  in  this  study. 

2.1  Equipment  Setup 

A  behavioral  science  engineering  laboratories  70  kpa  triaxial  test  cell  was 
used  to  apply  the  lateral  confining  pressure.  The  test  cell  had  a  pressure 
transducer  that  was  used  to  monitor  fluid  pressure  inside  the  cell.  The  axial 
load  was  applied  with  a  universal  testing  machine. 

Measurements  Group,  Inc.,  EA-06-500BL-350  strain  gauges  were  used  in 
56%  of  the  tests  (Wang  and  Ruthland,  1995).  Only  load-displacement  data 
were  recorded  in  the  remaining  tests.  Axial  load  and  displacement  data  were 
recorded  on  an  X-Y  plotter.  Some  of  the  test  data  were  digitized  and  stored  on 
an  IBM-compatible  personal  computer  using  an  analog  to  digital  converter. 
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2.2  Specimen  Preparation 

All  triaxial  specimens  were  prepared  with  a  water/cement  ratio  between  .42 
and  .47.  The  cement/aggregate  ratio  was  .25  by  weight.  Superplasticizer,  Friz- 
pak  Supercizer  6,  was  added  to  increase  the  workability  of  the  mixes  .3%  by 
weight  of  cement. 

Seventy-five  specimens  measuring  5  cm  in  diameter  by  10  cm  in  height 
were  cast  in  three  batches  of  25.  Twelve  specimens  5  cm  in  diameter  by  20  cm 
in  height  were  cast  in  a  separate  batch. 

Cap  blocks  were  prepared  to  help  relieve  the  friction  between  the  platens 
and  the  specimens.  The  cap  blocks  were  5  cm  in  diameter  and  1.3  to  2.5  cm 
thick.  The  mix  design  for  the  cap  blocks  was  similar  to  that  of  the  specimens 
except  that  the  cement/aggregate  ratio  was  increased  to  .4,  the  water/cement 
ratio  was  decreased  to  .35,  and  the  amount  of  superplasticizer  was  increased 
to  .5%.  The  goal  was  to  produce  a  material  with  similar  deformation 
properties  but  with  a  higher  strength.  The  modulus  and  lateral  extension  ratio 
of  the  two  materials  were  similar  (within  5%),  while  the  uniaxial  strength  of 
the  cap  blocks  was  50%  larger  than  the  strength  of  the  specimens. 

All  specimen  and  cap  blocks  were  polished  on  a  lap  table  containing  a 
coarse  grit  (100)  lap  pad.  The  cylinders  and  cap  blocks  contained  air  voids  on 
the  surfaces  that  would  tear  the  triaxial  membranes.  In  order  to  protect  the 
membranes  the  surface  voids  were  filled  with  clay.  Triaxial  membranes  were 
placed  on  the  specimens. 

Several  methods  of  relieving  the  friction  between  the  specimen  and  the 
platens  were  tested.  Based  on  some  preparing  tests,  cap  blocks  alone,  cap 
blocks  with  grease,  and  cap  blocks  with  polyethylene  sheets,  were  used  to 
relieve  the  friction.  All  the  tests  were  performed  at  0,  1.7,  3.5,  7.0,  14.0,  28.0, 
42.0,  and  56.0  Mpa,  and  they  were  conducted  using  a  standard  triaxial  loading 
where  the  displacement  was  monotonically  increased  up  to  and  past  the  peak 
load. 

3  Mesostructure  Simulation 

In  meso-  or  micro-scopic  studies,  it  is  necessary  to  generate  a  realistic  3- 
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dimensional  multiphase  system  consisting  of  cement  paste  and  coarse 
aggregate,  especially  for  concrete  (Desiderio,  1996).  This  is  because 
mechanical  behaviors  are  determined  by  extremes  in  the  local  stresses,  which 
have  a  close  relationship  with  the  distribution  of  cement  pastes  and  coarse 
aggregates,  as  well  as  with  aggregate  sizes  in  the  meso-structure  simulation. 
Here,  the  Monte  Carlo  method  (Hassold,  1996)  is  used  to  simulate  the 
cylinder  specimen  measuring  5  cm  in  diameter  and  10  cm  in  height. 

Automatic  mesh  generation  produced  148,500  lattices  for  the  cylinder 
sample  (  Fig.  1).  Generally,  that  amount  of  lattices  is  necessary  to  reach  a 
properly  aggregated  size  of  the  given  concrete  specimen.  A  non-zero  index, 
which  corresponds  to  both  the  phase  and  orientation  of  the  aggregate,  is 
initialized  to  each  lattice  by  a  random  function.  In  a  2-phase  simulation  (a 
and  p ),  the  sign  of  S ,  indicates  the  phase  present  at  that  site,  while  the 

absolute  value  of  S ,  corresponds  to  the  orientation  of  the  grain  in  which  the 

site  is  embedded.  Sites  with  one  or  more  different  closest  neighbors  are 
interface  sites,  sites  with  only  like  nearest  neighbors  are  interior  sites.  Total 
system  energy  is  specified  by  assigning  a  positive  energy  to  interface  areas  (or 
cement  pastes)  and  zero  energy  to  the  interior  sites  (or  coarse  aggregates), 
which  can  be  calculated  via  the  Hamiltonian  equation  (Hockin,  1980) 

t  {  (J„+JM+2J^)[l-S(S,^,)]+(J„+J„-2Jrt) 

8  i=l  j=  1 

Sign(S ,  )Sign(S ,  )+(J„  -J„  )[Sign(S ,  )+Sign(S , )]  }  (1) 

In  equation  (1),  N  is  the  total  sites  of  system.  Z  is  the  number  of  neighbors 
of  ith  site.  8  is  the  Kronecker  delta  function.  ‘Sign’  is  the  sign  function 
defined  as  Sign(  S  t )  equals  1  if  S ,  is  positive,  and  Sign(  S  i )  equals  -1  if 
S is  negative.  Joa ,  J  pp  and  Jap  are  a  -a  ,  p  -  p ,  and  a  -  p 
interfacial  energies,  respectively. 

In  iteration,  the  kinetic  energy  of  the  aggregate  is  adjusted  by  the  Monte 
Carlo  technique.  First,  a  lattice  site  and  a  site  index  are  chosen  at  random.  The 
index  of  the  chosen  site  is  then  changed  to  the  new  index  if  the  corresponding 
total  system  energy,  H,  does  not  increase.  After  each  attempted  index  change, 
the  Monte  Carlo  step  increases,  and  a  simulated  mesostructure  is  yielded.  Fig. 
2  is  a  3-dimensional  simulated  mesostructure  plot  of  the  cylinder  specimen 
that  involves  148,500  lattices.  Different  types  of  aggregates  are  composed  of 
certain  numbers  of  lattices,  and  between  the  aggregates  are  the  interfaces  or 
cement  pastes. 
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Fig.l  Specimen  lattice  Plot 


Fig.2  Specimen  aggregate  distribution 


4  Model  Description 

Simulations  of  the  static  tests  were  performed  using  the  finite  element  analysis 
code,  which  was  developed  for  static  analysis  of  nonlinear  structural  systems. 
The  Drucker-Prager  criterion  and  the  low  tensile  strength  criterion  are  used. 
Based  on  comparison  of  the  experimental  and  numerical  results,  the  strength 
ratio  between  the  cement  paste  (including  interface)  and  coarse  aggregate  was 
randomly  distributed  between  0.1  to  0.4.  Table  1  is  the  parameter  used  in  the 
analysis. 


Table  1  Parameters  for  coarse  aggregate* 


Young’s 

Modulus 

Poisson’s  Ratio 

Maxi.  Tensile 
Stress 

Maxi.  Comp¬ 
ressive  Stress 

Cohesion 

17  MPa 

0.167 

20  MPa 

100  Mpa 

22.4  MPa 

*Strength  ratio  between  cement  paste  (including  interface)  and  aggregate:  0.1  ~0.4 
randomly. 


An  8-node  isoparameter  3-dimensional  element  was  used.  In  this  case,  there 
were  148,500  lattices  in  the  cylinder  model  and  each  of  them  possessed 
different  mechanical  properties,  such  as  ultimate  strength  and  modulus  of 
elasticity,  which  need  to  be  incorporated  in  numerical  analyses.  However,  it  is 
very  difficult  or  almost  impossible  to  treat  that  many  elements  in  finite 
element  analysis  due  to  the  low  efficiency  and  floating  point  errors  of  the 


calculation.  The  situation  becomes  even  more  difficult  with  nonlinear 
analysis. 

Modified  Gauss  point  method  has  been  successfully  used  to  reduce  the 
computer  processing  time  while  still  guaranteeing  accuracy  (  Wang  and  Chen 
,1997).  The  basic  idea  of  the  modified  Gauss  point  method  used  here  is  that 
each  Gauss  point  in  the  numerical  integration  corresponds  to  one  lattice  and, 
therefore,  one  element  will  involve  a  specific  number  of  lattices,  say  8  or  27 
in  a  3-dimensional  problem.  The  Gauss  point  coordinates  and  the 
corresponding  weight  coefficients  need  to  be  modified  because  of  the  evenly 
distributed  lattices.  Table  2  lists  the  values  of  Gauss  point  coordinates,  ,  and 
the  weighted  coefficient,  Hj  for  the  different  number  of  Gauss  points,  n, 
evenly  distributed  between  -1  and  +1  in  case  of  1 -dimensional. 


Table  2  Values  of  corresponding  to  various  n 


n 

Hi 

6 

Precise  for  mth  Order 
Polynomial 

3 

2/3 

±0.707107  0.0 

m=3 

m 

±0.794654  ±0.187593 

m=4  or  5 

m 

2/5 

±0.832497  0.0  ±0.374541 

m=5 

6 

2/6 

±0.866246  ±0.422522  ±0.374541 

m=5  or  6 

-0.005  0  0.005  0.01  0.015  0.02  0.025 

DISPLACEMENT 


Fig. 3(a)  Typical  displacement-stress  relationship  of  experimental  data  (Units  of 

displacement  is  in  cm) 
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Displacement  vs  Stress 


Displacement 

Fig.3(b)  Typical  displacement-stress  relationship  Numerical  simulation  (Units  of 

displacement  is  in  cm) 


In  this  paper,  27  modified  Gauss  points  or  lattices  were  involved  in  one  8- 
node  element.  Therefore,  the  cylinder  specimen  had  5,500  elements,  which 
stand  for  148,500  lattices.  Several  typical  examples  (Wang  and  Chen,  1997) 
show  that  the  modified  Gauss  point  method  not  only  reduces  the  calculating 
time  greatly  but  also  guarantees  the  accuracy  of  the  computation. 


5  Results  and  Discussion 

Results  of  the  experimental  measurements  and  numerical  analyses  are  shown 
in  Fig.  3  through  5.  Very  good  agreement  is  found  between  displacement  and 
stress  (Fig.  3(a,  b)).  Fig.  4(a)  shows  the  failure  specimens  by  test,  while  Fig. 
4(b)  shows  the  damaged  cylinders  by  numerical  simulation  at  various 
confining  pressures.  Fig.  5  represents  comparison  of  the  angles  of  the  failure 
planes  by  experimental  and  numerical  analyses.  It  shows  that  a  distinct  failure 
surface  or  surfaces  change  orientation  from  vertical  to  angled  as  confinement 
increases.  Based  on  this  research,  there  are  several  points  that  need  to  be 
emphasized.  First,  confining  pressure  due  to  friction  between  platen  and 
specimen  has  a  predominating  effect  on  the  failure  angle  of  a  specimen  under 
triaxial  loading.  It  is  essential  to  relieve  the  friction  between  the  platens  and 
the  specimens.  Second,  to  predict  the  failure  mechanism  accurately,  randomly 
distributed  cement  paste  and  coarse  aggregate  is  necessary  in  the  simulation 
of  3-dimensional  concrete  specimen  under  triaxial  loading. 


29 


P  =  0.0  Mpa 


P  =  3.5  Mpa 


P  =  14.0  Mpa 

Fig.4(b)  Damaged  specimen  by 
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ANGLE  (d«gre«) 


Fig.4(a).  Damaged  specimen  by  experiments 


ANGLE  OP  THE  FAILURE  PLANE  VERSUS  CONFINING  PRESSURE 
The  angles  an  referenced  to  toe  direct**  of  maximum  compnssive  stress,  X j 
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Angle  versus  Confining  Pressure  by  FEM 


Angle  versus  confining  pressure  by  test(a)  and  numerical  simulation(b) 


6  Conclusion 


This  study  presented  a  procedure  to  predict  the  failure  mechanisms  of 
concrete  under  triaxial  loading.  Both  experimental  and  numerical  studies  showed 
changes  in  the  orientation  of  the  failure  plane  as  confining  pressure  increased  for 
concrete  cylinder  specimens.  For  realistic  triaxial  loading  in  experiments,  it  is 
essential  to  relieve  the  friction  between  the  platen  and  specimen.  Based  on 
comparison  of  experimental  and  numerical  results,  the  strength  ratio  between 
cement  paste  (including  interface)  and  coarse  aggregate  is  randomly  distributed 
between  0.1  to  0.4.  It  is  strongly  suggested  that  a  realistic  3-dimensional  meso- 
structure  simulation  be  used  in  a  numerical  model.  The  presented  modified 
Gauss  point  method  is  a  key  step  to  guarantee  the  efficiency  and  accuracy  of  the 
simulation. 
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One  of  the  major  technical  interests  for  quasi-brittle  materials  is  the  failure 
mechanism  under  uniaxial  or  triaxial  loading,  or  under  some  other  standard 
tests  such  as  beam  bending  specimen  and  compact  specimen.  Recent  research 
has  led  to  development  of  formal  procedures  for  predicting  damage  behavior  of 
quasi-brittle  material  using  computational  simulation  in  the  tests  mentioned 
above.  These  procedures  have  subsequently  been  programmed  into  a  computer 
module  and  thus  have  helped  the  researcher  have  a  general  direction  of  what 
the  corresponding  experiment  will  be.  The  objective  of  this  paper  is  to  present 
and  describe  results  obtained  using  meso-  or  micro-mechanical  simulation  of 
quasi-brittle  materials  such  as  concrete  and  ceramic.  Typical  results  show  that 
the  presented  method  provides  a  less  expensive  alternate  to  understanding  the 
failure  mechanisms  of  quasi-brittle  materials. 

Keywords:  failure  mechanism;  quasi-brittie;  structure  simulation; 
modified  Gauss  point 


Introduction 

While  the  age  of  quasi-brittle  materials  including  concrete  and  ceramic  is  upon  us, 
the  ability  to  accurately  predict  their  performance  mechanically  has  not  yet  arrived. 
In  view  of  the  difficulty  in  obtaining  the  failure  mechanisms  of  these  materials,  one 
must  often  rely  on  both  experimental  and  computational  simulation,  the  latter 
provides  a  much  less  expensive  alternative  to  the  former.1 

Concrete  and  ceramics,  being  typical  quasi-brittle  materials,  are  widely  used  in 
advanced  scientific  experiments  and  applications.  Computational  simulation  of 
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some  of  their  standard  tests  gives  the  researcher  an  easy  and  less  expensive  way  to 
master  their  failure  mechanisms.  Prior  to  numerical  analyses,  it  is  essential  to 
establish  a  realistic  3-dimensional  meso-  or  micro-structure  simulation.  Generally, 
hundreds  and  thousands  of  lattices  or  elements  need  to  be  included  in  the  simulated 
model.  This  is  extremely  important  because  mechanical  behaviors  are  determined 
by  extremes  in  the  local  stresses,  which  have  a  close  relationship  with  the 
distribution  of  grains  and  grain  boundaries  as  well  as  grain  size2.  Thus,  in  order  to 
simulate  the  ceramics  realistically,  both  in  a  geometrical  and  mechanical  sense, 
enough  simulated  elements  or  lattices  are  essential. 

In  this  paper,  the  Monte  Carlo  method  is  used  to  simulate  the  3-dimensional 
microstructure.  Because  numerical  analyses  of  a  meso-  or  microstructure  consume 
a  lot  of  computer  processing  time,  a  modified  Gauss  point  method  is  applied  to 
reduce  the  processing  time  while  still  guaranteeing  accuracy. 


Meso-  or  Micro-Structure  Simulation 


In  meso-  or  micro-scopic  studies  of  quasi-brittle  material  such  as  concrete,  it  is 
necessary  to  generate  a  realistic  3-dimensional  multiphase  system  consisting  of 
cement  paste  and  coarse  aggregate3.  As  an  example,  a  concrete  cylinder  specimen 
(Figure  1)  measuring  5  cm  in  diameter  and  10  cm  in  height  is  divided  into  148,500 
lattices.  Generally,  that  amount  of  lattices  is  necessary  to  reach  a  properly 
aggregated  size  of  the  given  concrete  specimen.  Afterward,  a  non-zero  index, 
which  corresponds  to  both  the  phase  and  orientation  of  the  aggregate,  is  initialized 
to  each  lattice  by  a  random  function.  In  a  2-phase  simulation  (a  and  (3 ),  the  sign 
of  Sj  indicates  the  phase  present  at  that  site,  while  the  absolute  value  of  Si 

corresponds  to  the  orientation  of  the  grain  in  which  the  site  is  embedded.  Sites 
with  one  or  more  different  closest  neighbors  are  interface  sites,  sites  with  only  like 
nearest  neighbors  are  interior  sites.  Total  system  energy  is  specified  by  assigning  a 
positive  energy  to  interface  areas  (or  cement  pastes)  and  zero  energy  to  the  interior 
sites  (or  coarse  aggregates),  which  can  be  calculated  via  the  Hamiltonian  equation 

4 

H=jl  £  {  (Jm+Jw+2J„p)[l-8(S,,SJ)]+(J„+Jpr2J,#) 

°  i=l  j=  1 

Sign(S  )Sign(S  j  )+(Jaa  -J  pp  )[Sign(S )+Sign(S  j )]  }  (1) 
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Figure  2  Specimen  Aggregate 


Figure  1  Specimen  Lattice  Plot 
Distribution 


In  equation  (1),  N  is  the  total  sites  of  the  system.  Z  is  the  number  of  neighbors 
of  the  ith  site.  8  is  the  Kronecker  delta  function.  ‘Sign’  is  the  sign  function  defined 
as  Sign(S  . ).  It  equals  1  if  S  ( is  positive,  and  Sign(S  i )  equals  -1  if  S ;  is  negative. 
Jaa ,  Jpp  and  Jap  are  a  -  a  ,  P  -  p ,  and  a  -  P  interfacial  energies,  respectively. 

In  iteration,  the  kinetic  energy  of  the  aggregate  is  adjusted  by  the  Monte  Carlo 
technique.  First,  a  lattice  site  and  a  site  index  are  chosen  at  random.  The  index  of 
the  chosen  site  is  then  changed  to  the  new  index  if  the  corresponding  total  system 
energy,  H,  does  not  increase.  After  each  attempted  index  change,  the  Monte  Carlo 
step  increases,  and  a  simulated  mesostructure  is  yielded.  Figure  2  is  a  3- 
dimensional  simulated  mesostructure  plot  of  the  cylinder  specimen  that  involves 
148,500  lattices.  Different  types  of  aggregates  are  composed  of  certain  numbers  of 
lattices,  and  between  the  aggregates  are  the  interfaces  or  cement  pastes. 

The  same  procedures  work  for  ceramic  samples  that  include  grains  and  their 
interfaces. 


Model  Description 

Simulations  of  the  static  tests  were  performed  using  the  finite  element  analysis 
code,  which  was  developed  for  static  analysis  of  nonlinear  structural  systems.  The 
Drucker-Prager  criterion  and  the  low  tensile  strength  criterion  are  applied. 

An  8-node  isoparameter  3-dimensional  element  is  used  in  this  code.  Generally, 
the  simulated  model  includes  hundreds  and  thousands  of  lattices,  each  of  them 
possessing  different  mechanical  properties,  such  as  ultimate  strength  and  modulus 
of  elasticity,  which  need  to  be  incorporated  in  numerical  analyses.  However,  it  is 
very  difficult  or  almost  impossible  to  treat  that  many  elements  in  finite  element 
analysis  due  to  the  low  efficiency  and  floating  point  errors  of  the  calculation.  The 
situation  becomes  even  more  difficult  with  nonlinear  analysis. 

Modified  Gauss  point  method  has  been  successfully  used  to  reduce  the 
computer  processing  time  while  still  guaranteeing  accuracy5.  The  basic  idea  of  the 
modified  Gauss  point  method  used  here  is  that  each  Gauss  point  in  the  numerical 
integration  corresponds  to  one  lattice  and,  therefore,  one  element  involves  a 
specific  number  of  lattices,  say  8  or  27  in  a  3-dimensional  problem.  The  Gauss 
point  coordinates  and  the  corresponding  weight  coefficients  need  to  be  modified 
because  of  the  evenly  distributed  lattices.  Table  1  lists  the  values  of  Gauss  point 
coordinates,  ,  and  the  weighted  coefficient.  Hi,  for  the  different  number  of  Gauss 

points,  n,  evenly  distributed  between  -1  and  +1  in  the  case  of  1-dimensional 
problem. 


Table  1  Values  of  t)j  Corresponding  to  Various  n 


n 

Hi 

& 

Precise  for  mth  Order 
Polynomial 

3 

2/3 

±0.707107  0.0 

m=3 

4 

2/4 

±0.794654  ±0.187593 

m=4  or  5 

5 

2/5 

±0.832497  0.0  ±0.374541 

m=5 

6 

2/6 

±0.866246  ±0.422522  ±0.374541 

m=5  or  6 

In  this  paper,  27  modified  Gauss  points  or  lattices  are  involved  in  one  8-node 
element.  As  an  example,  the  cylinder  specimen  shown  in  Figure  1  has  5,500 
elements  which  stand  for  148,500  lattices.  Numerical  examples  show  that  the 
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modified  Gauss  point  method  not  only  reduces  the  calculating  time  greatly  but  also 
guarantees  the  accuracy  of  the  computation. 


Numerical  Examples  and  Discussions 


1  Four  Points  Beam  Bending 


The  chosen  numerical  example  is  the  so-called  four  points  beam  bending  problem 
(Figure  3).  The  beam  has  a  cross  section  of  3  cm  (height)  x  1  cm  (width),  and  a 
span  of  18  cm.  The  Young’s  modulus  is  200  GPa,  the  Poisson  ratio  is  0.2,  and  the 
cohesive  strength  is  17  Mpa.  The  material  of  the  beam  is  treated  as  a  realistic 
elasto-plastic.  Two  models  are  taken  into  consideration: 
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Figure  3  4-Points  Beam  Bending  Specimen 


(1)  Common  Gauss  Point  Method 

3x3x3  Gauss  points  are  used  in  one  element  to  carry  the  numerical  integration. 
A  total  of  216  elements  are  included  in  this  model. 


(2)  Modified  Gauss  Point  Method 

3x3x3  lattices  are  included  in  each  element.  A  total  of  2 16  x  27,  or  5832  lattices 
are  included  in  the  second  model. 

The  above  two  models  produce  the  same  numerical  results.  The  yield  zoom  is 
shown  in  Figure  3  that  occurs  at  P  =  2126  N.  Moreover,  the  two  models  consume 
almost  the  same  amount  of  run  time  although  the  second  one  can  simulate  more 
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lattices  than  the  first.  Table  2  shows  the  comparison  between  the  theoretical  and 
numerical  results. 


Table  2:  Comparison  of  Horizontal  Stresses  at  h/4  in  the  Yield  Zoom 


Theoretical 

Numerical 

Horizontal  stress 

29.45  MPa 

29.31  MPa 

This  example  suggests  that  although  approximately  the  same  amount  of  run  time 
is  consumed  because  they  include  the  same  number  of  integral  Gauss  points,  the 
one  using  the  modified  Gauss  point  method  could  consider  much  more  lattices  and 
produce  satisfactory  result  at  the  same  time. 


2  Cylinder  Specimen  under  Triaxial  Loading 


This  example  describes  the  failure  mechanisms  of  a  concrete  specimen  under 
triaxial  loading.  The  failure  is  characterized  by  a  distinct  failure  surface  or  surfaces 
that  change  orientation  from  nearly  vertical  to  angled  as  confining  pressure 
increase.6  The  simulated  model  that  has  148,500  lattices  is  shown  in  Figure  2. 
After  applying  the  modified  Gauss  point  method,  the  specimen  involves  5,500 
elements. 

Uniaxial  displacement  with  a  5  pm  in  step  is  applied  on  the  top  of  the  specimen 
while  a  certain  amount  of  confining  pressure  is  exerted  on  the  surrounding  surface 
of  the  cylinder.  Based  on  the  results  of  the  experiment,6  the  friction  between  the 
platens  and  the  specimen  is  relieved  since  this  friction  severely  affects  the 
development  of  the  failure  mechanisms  as  confining  pressure  increases. 

Results  of  numerical  analyses  are  shown  in  Figure  4  through  6.  Very  good 
agreement  is  found  between  displacement  and  stress  (Figure  4(a,  b)).  Figure  5(a) 
shows  the  failure  specimens  by  test,  while  Figure  5(b)  shows  the  damaged 
cylinders  by  numerical  simulation  at  various  confining  pressures.  Figure  6 
represents  comparison  of  the  angles  of  the  failure  planes  by  experimental  and 
numerical  analyses.  It  shows  that  a  distinct  failure  surface  or  surfaces  change 
orientation  from  vertical  to  angled  as  confinement  increases.  Based  on  this 
research,  there  are  several  points  that  need  to  be  emphasized.  First,  confining 
pressure  due  to  friction  between  platen  and  specimen  has  a  predominating  effect 
on  the  failure  angle  of  a  specimen  under  triaxial  loading.  It  is  essential  to  relieve 
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the  friction  between  the  platens  and  the  specimens.  Second,  to  predict  the  failure 
mechanism  accurately,  randomly  distributed  cement  paste  and  coarse  aggregate  is 
necessary  in  the  simulation  of  3-dimensional  concrete  specimen  under  triaxial 
loading.  Finally,  to  achieve  the  best  results,  simple  and  reliable  nonlinear  criteria 
should  be  chosen.  For  simplicity,  criteria  with  easily  determined  parameters  are 
recommended. 


(a)  (b) 

Figure  4  Typical  Displacement-Stress  Relationship:  (a)  Experimental  Data  and 
(b)  Numerical  Simulation  (Units  of  Displacement  is  in  cm) 


Figure  5(a).  Damaged  Specimen  by  Experiments 


P  =  O.OMpa  P=1.7Mpa  P  =  3.5  Mpa 


P  =  7.0  Mpa  P  =  14.0  Mpa  P  =  28.0  Mpa 


Figure  5(b)  Damaged  Specimen  by  Numerical  Simulation 
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Angle  versus  Confining  Pressure  by  FEM 


Figure  6  Angle  versus  confining  pressure  by  test(a)  and  numerical  simulation(b) 


3  Compact  Tensile  Specimen 


A  compact  tensile  specimen  was  analyzed  in  this  example.  The  configuration  of  the 
specimen  was  25mm( width)  x30mm(height)  xl.5mm(depth).  The  specimen  was 
loaded  by  a  load-point  displacement  that  is  separated  into  12  steps  with  2.5pm 
each.  A  modified  Gauss  point  method  is  applied  in  both  of  the  following  cases. 

(1)  The  Effect  of  Grain  Size 

Two  finite  element  models  are  put  into  consideration  with  total  lattice  numbers  of 
6701  and  15,182,  respectively.  Figures  7  and  8  show  the  grain  size  distributions 
with  800  and  1100  pm,  respectively.  Figure  9  is  the  corresponding  Peq  vs 
displacement  plot  that  reveals  that  the  smaller  the  grain  size,  the  smaller  the 
compact  tensile  Peq.  This  conclusion  is  the  same  as  mentioned  in  reference  7. 

(2)  The  Effect  of  Boundary  Conditions 

Fixed  and  free  boundary  conditions  are  applied  to  the  above  specimen(Figure  10). 
This  model  has  6701  lattices  and  corresponding  grain  size  is  1100  pm  as  shown  in 
Figure  7.  Figure  11  is  the  plot  of  equivalent  load  Peq  vs  applied  known 
displacement,  where  Peq  is  calculated  by  the  equivalent  condition.  Figures  12 
shows  the  comparison  of  the  final  damaged  specimen  with  two  different  boundary 
conditions.  According  to  the  profiles  of  each  step  in  simulation,  the  cracking  starts 
from  the  crack  tip  and  develops  as  the  applying  tensile  force  increases.  Several 
single  cracks  occur  during  the  simulation.  Some  of  them  develop  while  others  stop 
cracking  or  even  go  back  to  closure.  The  final  failure  does  not  take  place  until  one 
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or  two  main  cracks  reach  the  critical  cracking  length.  Fig.  13-15  are  the 
experimental  profiles  of  the  cracking  development  that  match  the  numerical 
simulation  very  well. 


Grain  Size  film)  Crain  Sit' (»ml 


Figure  7  Grain  Size  Distribution  Figure  8  Grain  Size  Distribution 


0.0  5.0  10.0  15.0  20.0  25.0  30.0 

Displacement,  <Mm) 

Figure  9  Effort  s  of  Grain  Size 
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Free  Boundary  Condition  Fixed  Boundary  Condition 

Figure  10  Configuration  of  Specimen 


Figure  1 1  Effect  of  Boundary  condition  (Grain  Size:  1 100  |im) 
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(a)  (b) 

Figure  12  Damaged  Specimen  (a)  Fixed  Boundary  (b)  Free  Boundary  Condition 
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Fig.  13  Microscope  Profile  of  Cracking  Specimen  (Step  1 ) 
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Fig.  14  Microscope  Profile  of  Cracking  Specimen  (Step  2  ) 
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Fig.  15  Microscope  Profile  of  Cracking  Specimen  (Step  3  ) 


Conclusion 


This  study  presented  a  procedure  to  predict  the  failure  mechanisms  of  quasi-brittle 
materials  under  triaxial  loading  and  some  other  standard  tests.  Prior  to  numerical 
analyses,  it  is  essential  to  have  a  realistic  3-dimensional  meso-  or  micro-structure 
simulation.  The  presented  modified  Gauss  point  method  is  another  key  step  to 
guarantee  the  efficiency  and  accuracy  of  the  simulation.  The  mechanical  properties 
and  failure  mechanisms  predicted  by  the  above  method  are  in  good  agreement  with 
the  previous  experimental  results.  Thus,  the  method  provides  a  less  expensive 
alternative  in  understanding  the  failure  mechanisms  of  quasi-brittle  materials. 
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